function solve_counter = solve_counter(fp,param)
alpha = param(1:9);
TFP = param(10:11) ;  

init_draws      = cell(4,1);
peers_forward   = cell(4,1);
peers_backward  = cell(4,1);
PS_forward      = cell(4,1);

h = 1    ;

    for j = [fp.origin_neigh fp.receiving_neigh]

        init_draws{j}     = fp.schools_distrib(j,1) + fp.schools_distrib(j,2).*(fp.init_draws{j}(1:fp.schools_distrib(j,3),1,h) - mean(fp.init_draws{j}(1:fp.schools_distrib(j,3),1,h)))./std(fp.init_draws{j}(1:fp.schools_distrib(j,3),1,h));    

        %Shocks;
        if j==fp.origin_neigh

            peers_forward{j}  = fp.peers_forward{j}(1:fp.schools_distrib(j,3),1:fp.schools_distrib(j,3),h,:);
            peers_backward{j} = fp.peers_backward{j}(:,1:fp.schools_distrib(j,3), :,:);                           
            PS_forward{j}     = fp.PS_forward{j}(1:fp.schools_distrib(j,3),h,:);    
    
        else
            
            peers_forward{j}  = fp.peers_forward{j}(1:(fp.schools_distrib(j,3)+fp.n_kids_moved),1:(fp.schools_distrib(j,3)+fp.n_kids_moved),h,:);
            peers_backward{j} = fp.peers_backward{j}(:,1:(fp.schools_distrib(j,3)+fp.n_kids_moved), :,:);    
            PS_forward{j}     = fp.PS_forward{j}(1:(fp.schools_distrib(j,3)+fp.n_kids_moved),h,:);                
            
            
        end
    end
        
    %Ordering children to be moved from the sending neighborhoods;
    temp_init_draws = sort(init_draws{fp.origin_neigh}(:,1), 1, 'ascend');
    
    
    if fp.do_counter_type==1
       
    if fp.n_kids_moved ==1    
    %Moving the Median Child;    
    moved_children_skills = temp_init_draws( round( fp.schools_distrib(fp.origin_neigh,3)/2 ) , 1);      
    end
    
    if fp.n_kids_moved==0
    moved_children_skills  = [];
    end
    end
    
    
    if fp.do_counter_type==2 || fp.do_counter_type==3 
    
    %Gradually moving children;
    moved_children_skills = temp_init_draws(fp.shock_child_moved(:,h) , 1 ) ;
        
    end
            
    
    %%%%%
    %Gathering/Ordering Data for Sent and Receiving Children in the Counterfactual Economy;
    %%%%%
    
    %%%%%%%%%%%%%%%%;
    %Allow for families from receiving neighborhood to leave;
    %%%%%%%%%%%%%%%%;

    %Calculating change in skills implied by the policy;
    temp_initial_distrib = [ init_draws{fp.receiving_neigh}; moved_children_skills ];
    temp_initial_mean = mean(temp_initial_distrib);
    temp_change_mean = (temp_initial_mean -   fp.schools_distrib(fp.receiving_neigh,1));


    prob_flight_neighborhood = zeros(fp.schools_distrib(4,3),1);

    for x = 1:1:10
        
        prob_flight_neighborhood(init_draws{fp.receiving_neigh}>= fp.deciles_initial_skills(x) &   init_draws{fp.receiving_neigh}<fp.deciles_initial_skills(x+1),1) = 100.*temp_change_mean.*fp.elasticity_moving_by_kid_skills(x);        

    end

    %Simulate families that decide to leave the neighborhood because of the
    %policy;
    flight_neighborhood(:,1) = (fp.shocks_leaving_neighborhood(:,h)<=prob_flight_neighborhood);
    final_init_draw_receiving = init_draws{fp.receiving_neigh};
    final_init_draw_receiving(flight_neighborhood==1)=NaN; 

    %Initial Draws (Keep the same shock for every child Baseline vs
    %Counterfactual);
    fp_parfor.init_draws = [ final_init_draw_receiving ; moved_children_skills ];
    fp_parfor.peers_forward  = peers_forward{fp.receiving_neigh};
    fp_parfor.peers_backward = peers_backward{fp.receiving_neigh};    
    fp_parfor.PS_forward     = PS_forward{fp.receiving_neigh};    
    
    %Indicator if the child is moved or receiving;
    fp_parfor.moved_children = [ zeros(length(init_draws{fp.receiving_neigh}),1) ; ones( length(moved_children_skills) , 1)];
    
    fp_parfor.n_simulation = length(fp_parfor.init_draws);
    
    epsilon = 5;
    toll= 0.05;

    record_epsilon=[3 epsilon];
    
    epsilon_skills = 5;
    
    i = 1;
    j = fp.receiving_neigh;


if fp.do_counter_type < 3    
    
while (epsilon > toll  && epsilon_skills>toll)  %Fix point on the policy functions;
    
fp_parfor.i = i; %iteration   
fp_parfor.j = j; %neighborhood

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;    
%Simulate Forward;    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

if i == 1

    %first iteration: set the initial guess about the distribution of
    %skills;    
        
    fp_parfor.distrib_PS = zeros( size( fp_parfor.init_draws , 1 ) ,fp.periods-1);     

    fp_parfor.inv0{1} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    fp_parfor.inv0{2} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    fp_parfor.inv0{3} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    
    fp_parfor.inv1{1} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    fp_parfor.inv1{2} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    fp_parfor.inv1{3} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);    
    
    distrib_skills1 = exp(fp_parfor.init_draws(:,1));            
    distrib_skills2 = technology(1,alpha, distrib_skills1 , fp.Max_Time, fp_parfor.distrib_PS(:,1) ,  distrib_skills1 , TFP) ;
    distrib_skills3 = technology(2,alpha, distrib_skills2 , fp.Max_Time, fp_parfor.distrib_PS(:,2) ,  distrib_skills2 , TFP)  ;
    distrib_skills4 = technology(3,alpha, distrib_skills3 , fp.Max_Time, fp_parfor.distrib_PS(:,3) ,  distrib_skills3 , TFP)  ;
              
    fp_parfor.distrib_skills(:,1) = distrib_skills1;
    fp_parfor.distrib_skills(:,2) = distrib_skills2;
    fp_parfor.distrib_skills(:,3) = distrib_skills3;
    fp_parfor.distrib_skills(:,4) = distrib_skills4;    
    
    
    fp_parfor.grid{1} = prctile( distrib_skills1 , fp.percentiles_kid ) ;
    fp_parfor.grid_peers{1} = fp_parfor.grid{1} ;    
    fp_parfor.grid{2} = prctile( distrib_skills2 , fp.percentiles_kid ) ;
    fp_parfor.grid_peers{2} = fp_parfor.grid{2} ;              
    fp_parfor.grid{3} =  prctile(distrib_skills3 , fp.percentiles_kid );
    fp_parfor.grid_peers{3} = fp_parfor.grid{3};

else
    %solve forward to simulate the distribution of skills and investments;
    sf = solve_forward_counter(param,fp,fp_parfor,h);
    previous_iter_distrib_skills = fp_parfor.distrib_skills;    

    fp_parfor.distrib_skills = sf.distrib_skills;
    fp_parfor.distrib_PS = sf.distrib_PS;
    
    fp_parfor.grid{1} = linspace( prctile(exp(fp_parfor.init_draws(:,1)) , 5) , prctile(exp(fp_parfor.init_draws(:,1)), 95), length(fp.percentiles_kid) );
    fp_parfor.grid_peers{1} = fp_parfor.grid{1} ;
    fp_parfor.grid{2} = linspace( prctile( fp_parfor.distrib_skills(:,2) , 5) , prctile(fp_parfor.distrib_skills(:,2), 95 ), length(fp.percentiles_kid) );
    fp_parfor.grid_peers{2} = fp_parfor.grid{2} ;              
    fp_parfor.grid{3} = linspace( prctile( fp_parfor.distrib_skills(:,3), 5 ) , prctile(fp_parfor.distrib_skills(:,3), 95 ), length(fp.percentiles_kid) );
    fp_parfor.grid_peers{3} = fp_parfor.grid{3} ;    
 
   fp_parfor.inv0{1} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
   fp_parfor.inv0{2} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
   fp_parfor.inv0{3} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    
   fp_parfor.inv1{1} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
   fp_parfor.inv1{2} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
   fp_parfor.inv1{3} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);      
    
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Solve Backward;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

sb= solve_backward_counter(param, fp,fp_parfor) ;      

fp_parfor.Hc_tp1_guess_prime = sb.Hc_tp1_guess_prime;

if i > 1 
eps1=zeros(fp.periods-1,1);
eps2=zeros(fp.periods-2,1);
eps_skills=zeros(fp.periods-2,1);


for t = 1:fp.periods-1
            
  eps1(t) = max(abs(fp_parfor.estim_policy0{t}(:) - sb.estim_policy0{t}(:)));
  
  if t==3

  end
  
    if t < fp.periods-1
        
        eps2(t) = max(abs(fp_parfor.estim_policy1{t}(:) - sb.estim_policy1{t}(:)));
 
    end
    
   if t > 1  
   deciles_skills_previous_iter = quantile(log(previous_iter_distrib_skills(:,t)),[0.05:0.1:0.95]);
   deciles_skills_current_iter = quantile(log(sf.distrib_skills(:,t)),[0.05:0.1:0.95]);          
   
   eps_skills(t) =  max(abs(deciles_skills_previous_iter-deciles_skills_current_iter));
          
   end    
    
    
end

epsilon = max([eps1;eps2]);
epsilon_skills = max(eps_skills);
%keep track of epsilon records;
ind_rec = mod(i-1,2) + 1;
record_epsilon(ind_rec) = epsilon;


end

fp_parfor.estim_policy0 = sb.estim_policy0;
fp_parfor.estim_policy1 = sb.estim_policy1;

fp_parfor.estim_value0 = sb.estim_value0;
fp_parfor.estim_value1 = sb.estim_value1;


fprintf('%d ', i)

i = i + 1;
%record_epsilon
end

fprintf('Solved Neighborhood %d \n ', j)   

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Simulate forward once the fixed point is achieved and create the simulated 
%dataset;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
fprintf('Simulating Neighborhood %d:\n ', j )

for h = 1 : fp.n_rip_sim

if h ==fp.n_rip_sim
fprintf('Simulated Neighborhood %d \n ', j)
end

    for j = [fp.origin_neigh fp.receiving_neigh]

        init_draws{j}     = fp.schools_distrib(j,1) + fp.schools_distrib(j,2).*(fp.init_draws{j}(1:fp.schools_distrib(j,3),1,h) - mean(fp.init_draws{j}(1:fp.schools_distrib(j,3),1,h)))./std(fp.init_draws{j}(1:fp.schools_distrib(j,3),1,h));    

        %Shocks;
        if j==fp.origin_neigh

            peers_forward{j}  = fp.peers_forward{j}(1:fp.schools_distrib(j,3),1:fp.schools_distrib(j,3),h,:);
            PS_forward{j}     = fp.PS_forward{j}(1:fp.schools_distrib(j,3),h,:);    
    
        else
            
            peers_forward{j}  = fp.peers_forward{j}(1:(fp.schools_distrib(j,3)+fp.n_kids_moved),1:(fp.schools_distrib(j,3)+fp.n_kids_moved),h,:);
            PS_forward{j}     = fp.PS_forward{j}(1:(fp.schools_distrib(j,3)+fp.n_kids_moved),h,:);                
            
            
        end
    end
        
    
    temp_init_draws = sort(init_draws{fp.origin_neigh}(:,1), 1, 'ascend');
    
    
    if fp.do_counter_type==1
       
    if fp.n_kids_moved ==1    
    %Moving the Median Child;    
    moved_children_skills = temp_init_draws( round( fp.schools_distrib(fp.origin_neigh,3)/2 ) , 1);      
    end
    
    if fp.n_kids_moved==0
    moved_children_skills  = [];
    end
    end
    
    
    if fp.do_counter_type==2 || fp.do_counter_type==3 
    
    %Gradually moving children;
    moved_children_skills = temp_init_draws(fp.shock_child_moved(:,h) , 1 ) ;
        
    
    if fp.do_counter_type==3

    fp_parfor.j = j;
    
    data_origin_initial = fp.simul_data_baseline{fp.origin_neigh}( fp.simul_data_baseline{fp.origin_neigh}( : ,fp.ind_data.child_age,h)==9 ,:,h);
    data_origin_10      = fp.simul_data_baseline{fp.origin_neigh}( fp.simul_data_baseline{fp.origin_neigh}( : ,fp.ind_data.child_age,h)==10 ,:,h);
    data_origin_11      = fp.simul_data_baseline{fp.origin_neigh}( fp.simul_data_baseline{fp.origin_neigh}( : ,fp.ind_data.child_age,h)==11 ,:,h);
    temp_data           = sortrows( [data_origin_initial(:,[fp.ind_data.child_C   fp.ind_data.inv    fp.ind_data.ParStyle])  data_origin_10(:,[ fp.ind_data.inv    fp.ind_data.ParStyle]) data_origin_11(:,[ fp.ind_data.inv    fp.ind_data.ParStyle])   ] , 1) ;
    
    data_receiving_inv  =[  fp.simul_data_baseline{fp.receiving_neigh}( fp.simul_data_baseline{fp.receiving_neigh}( : ,fp.ind_data.child_age,h)==9  ,fp.ind_data.inv,h), ...
                            fp.simul_data_baseline{fp.receiving_neigh}( fp.simul_data_baseline{fp.receiving_neigh}( : ,fp.ind_data.child_age,h)==10 ,fp.ind_data.inv,h), ...
                            fp.simul_data_baseline{fp.receiving_neigh}( fp.simul_data_baseline{fp.receiving_neigh}( : ,fp.ind_data.child_age,h)==11 ,fp.ind_data.inv,h) ];
                        
    data_receiving_PS   =[  fp.simul_data_baseline{fp.receiving_neigh}( fp.simul_data_baseline{fp.receiving_neigh}( : ,fp.ind_data.child_age,h)==9  ,fp.ind_data.ParStyle,h), ...
                            fp.simul_data_baseline{fp.receiving_neigh}( fp.simul_data_baseline{fp.receiving_neigh}( : ,fp.ind_data.child_age,h)==10 ,fp.ind_data.ParStyle,h), ...
                            fp.simul_data_baseline{fp.receiving_neigh}( fp.simul_data_baseline{fp.receiving_neigh}( : ,fp.ind_data.child_age,h)==11 ,fp.ind_data.ParStyle,h) ];
    
    fp_parfor.fixed_inv = [ data_receiving_inv ; temp_data( fp.shock_child_moved(:,h) , [ 2 4 6 ] )]; 
    fp_parfor.fixed_PS  = [ data_receiving_PS  ; temp_data( fp.shock_child_moved(:,h) , [ 3 5 7 ] )];    
    
    end
    
    
    end
    
    id_moved_children(:,h) = ismember( init_draws{fp.origin_neigh}(:,1) , moved_children_skills);    
    
    
    %%%%%
    %Gathering/Ordering Data for Sent and Receiving Children in the Counterfactual Economy;
    %%%%%
    

    temp_initial_distrib = [ init_draws{fp.receiving_neigh}; moved_children_skills ];

    temp_initial_mean = mean(temp_initial_distrib);

    temp_change_mean = (temp_initial_mean -   fp.schools_distrib(fp.receiving_neigh,1));


    prob_flight_neighborhood = zeros(fp.schools_distrib(4,3),1);

    for x = 1:1:10
        
        prob_flight_neighborhood(init_draws{fp.receiving_neigh}>= fp.deciles_initial_skills(x) &   init_draws{fp.receiving_neigh}<fp.deciles_initial_skills(x+1),1) = 100.*temp_change_mean.*fp.elasticity_moving_by_kid_skills(x);        

    end

    flight_neighborhood(:,1) = (fp.shocks_leaving_neighborhood(:,h)<=prob_flight_neighborhood);

    final_init_draw_receiving = init_draws{fp.receiving_neigh};
    final_init_draw_receiving(flight_neighborhood==1)=NaN; 


    id_flights(:,h) = isnan(final_init_draw_receiving);


    %Initial Draws (Keep the same shock for every child Baseline vs
    %Counterfactual);
    fp_parfor.init_draws = [ final_init_draw_receiving ; moved_children_skills ];
    fp_parfor.peers_forward  = peers_forward{fp.receiving_neigh};
    fp_parfor.peers_backward = peers_backward{fp.receiving_neigh};    
    fp_parfor.PS_forward     = PS_forward{fp.receiving_neigh};    
    
    %Indicator if the child is moved or receiving;
    fp_parfor.moved_children = [ zeros(length(init_draws{fp.receiving_neigh}),1) ; ones( length(moved_children_skills) , 1)];
    
    fp_parfor.n_simulation = length(fp_parfor.init_draws);


sf = solve_forward_counter(param,fp,fp_parfor,h);



sim_data_counter(:,:,h) = sf.sim_data ; 
sim_data_network_counter(:,:,h) = sf.sim_data_network;
    
    

end


solve_counter.sim_data = sim_data_counter;
solve_counter.sim_data_network = sim_data_network_counter;
solve_counter.id_moved_children = id_moved_children;
solve_counter.id_flights        = id_flights;
